/*
* Copyright 1993-2010 NVIDIA Corporation.  All rights reserved.
*
* Please refer to the NVIDIA end user license agreement (EULA) associated
* with this source code for terms and conditions that govern your use of
* this software. Any use, reproduction, disclosure, or distribution of
* this software and related documentation outside the terms of the EULA
* is strictly prohibited.
*
*/

/* Template project which demonstrates the basics on how to setup a project 
* example application, doesn't use cutil library.
*/

#include <stdio.h>
#include <string.h>
#include <iostream>
#include <fstream>
#include <math.h>
#include <time.h>
#include "cuda.h"
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <cufft.h>
#include "fftw3.h"
#include "cudacomplex.h"


//#include <shrQATest.h>

using namespace std;

#define pi 3.141592653
#define NX 5
#define BATCH 1

struct Complex {
	float re;
	float im;
	__host__ __device__ Complex(float r, float i) {
		re = r;
		im = i;
	}
	__host__ __device__ Complex(const Complex& c) {
		re = c.re;
		im = c.im;
	}
	__host__ __device__ Complex& operator=(const Complex& c) {
		re = c.re;
		im = c.im;
		return *this;
	}
	__host__ __device__ Complex& operator*(const Complex& c) {
		float newR = re*c.re - im*c.im;
		float newI = re*c.im + im*c.re;
		re = newR;
		im = newI;
		return *this;
	}
	__host__ __device__ Complex& operator/(const Complex& c) {
		float div = c.re*c.re + c.im*c.im;
		float t_real = (re*c.re) + (im*c.im);
		t_real /= div;
		float t_imag = (im*c.re) - (re*c.im);
		t_imag /= div;
		im = t_imag;
		re = t_real;
		return *this;
	}
	__host__ __device__ Complex& operator+(const Complex& c) { 
		re = re + c.re;
		im = im + c.im;
		return *this;
	}
	__host__ __device__ Complex& operator-(const Complex& c) {
		re = re - c.re;
		im = im - c.im;
		return *this;
	}
	__host__ __device__ Complex& conjugate() {
		im = -im;
		return *this;
	}
};


int main(int argc, char **argv)
{

	if (argc != 0) {
		cout << "test" << endl;
		getchar();
	}

	int i;

	cufftHandle plan, planinv, plan2;
	cufftComplex data[NX*BATCH];
	cufftComplex *dev_data;

	/* source data creation */
	for(i=  0 ; i < NX*BATCH ; i++){
		data[i].x = (float) i;
		data[i].y = 0.0f;
	}
	for(i = 0 ; i < NX*BATCH ; i++){
		printf("data[%d] %f %f\n", i, data[i].x, data[i].y);
	}
	getchar();

	//printf("sizeof(cufftComplex) %d\n", sizeof(cufftComplex));

	cudaMalloc((void**)&dev_data, sizeof(cufftComplex)*NX*BATCH);
	/* transfer to GPU memory */
	cudaMemcpy(dev_data, data, sizeof(cufftComplex)*NX*BATCH, cudaMemcpyHostToDevice);

	/* Create a 1D FFT plan. */
	cufftPlan1d(&plan, NX, CUFFT_C2C, BATCH);
	/* Use the CUFFT plan to transform the signal in place. */
	cufftExecC2C(plan, dev_data, dev_data, CUFFT_FORWARD);

	/* transfer results from GPU memory */
	cudaMemcpy(data, dev_data, sizeof(cufftComplex)*NX*BATCH, cudaMemcpyDeviceToHost);

	/* Destroy the CUFFT plan. */
	cufftDestroy(plan);
	cudaFree(dev_data);

	for(i = 0 ; i < NX*BATCH ; i++){
		printf("data[%d] %f %f\n", i, data[i].x, data[i].y);
	}

	getchar();

	int N = 5;
	fftwf_complex *in, *out;
	fftwf_plan p,q;

	in = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex) * N);
	out = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex) * N);
	cudacomplex *in_Complex = (cudacomplex*) calloc(sizeof(cudacomplex), N);
	cudacomplex *out_Complex = (cudacomplex*) calloc(sizeof(cudacomplex), N);
	cudacomplex *out_cuda = (cudacomplex*) calloc(sizeof(cudacomplex), N);
	cudacomplex *out_cuda_inv = (cudacomplex*) calloc(sizeof(cudacomplex), N);
	cudacomplex *dev_in, *dev_out;
	cudaMalloc((void**)&dev_in, N * sizeof(cudacomplex));
	cudaMalloc((void**)&dev_out, N * sizeof(cudacomplex));


	for(i = 0 ; i < N ; i++){
		//printf("out[%d] %f %f\n", i, out[i][0], out[i][1]);
	}
	for(i=  0 ; i < N ; i++){
		in[i][0] = (float) i;
		//in_Complex[i] = make_cudacomplex((float) i, 0);
		in_Complex[i].real() = (float) i;
		//in[i][1] = 0.0f;
	}

	cudaMemcpy(dev_in, in_Complex, sizeof(cudacomplex)* N, cudaMemcpyHostToDevice);

	/* Create a 1D FFT plan. */
	int n[1];
	n[0] = N;
	int inembed[1];
	inembed[0] = N;
	int onembed[1];
	onembed[0] = 3;
	cufftPlanMany(&plan2, 1, n, inembed, 1, 1, onembed, 1, 3, CUFFT_C2C, BATCH );
	//cufftPlan1d(&plan2, NX, CUFFT_C2C, BATCH);
	//cufftPlan1d(&planinv, NX, CUFFT_C2C, BATCH);
	/* Use the CUFFT plan to transform the signal in place. */
	cufftExecC2C(plan2, (cufftComplex *) dev_in, (cufftComplex *) dev_out, CUFFT_FORWARD);
	cudaMemcpy(out_cuda, dev_out, sizeof(cudacomplex)* N, cudaMemcpyDeviceToHost);

	//cufftExecC2C(planinv, (cufftComplex *) dev_out, (cufftComplex *) dev_in, CUFFT_INVERSE);
	//cudaMemcpy(out_cuda_inv, dev_in, sizeof(cudacomplex)* N, cudaMemcpyDeviceToHost);
	

	/* Destroy the CUFFT plan. */
	cufftDestroy(plan2);
//	cufftDestroy(planinv);
	cudaFree(dev_in);
	cudaFree(dev_out);



	for(i = 0 ; i < N ; i++){
		printf("in_Complex[%d] %f %f\n", i, in_Complex[i].real(), in_Complex[i].imag());
	}		

	//p = fftwf_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
	p = fftwf_plan_dft_1d(N, (fftwf_complex*) in_Complex, (fftwf_complex*) out_Complex, FFTW_FORWARD, FFTW_ESTIMATE);
	q = fftwf_plan_dft_1d(N, (fftwf_complex*) out_Complex, (fftwf_complex*) in_Complex, FFTW_BACKWARD, FFTW_ESTIMATE);

	fftwf_execute(p); /* repeat as needed */
	fftwf_execute(q); /* repeat as needed */

	fftwf_destroy_plan(p);
	fftwf_destroy_plan(q);

	//cudacomplex mul;
	//mul = out_Complex[1] + out_Complex[2];
	//printf("real = %f\t imag = %f \n", mul.real(), mul.imag());

	//memcpy(out_Complex + 3, out_Complex, 2 * sizeof(cudacomplex));

	for(i = 0 ; i < N ; i++){
		printf("out_Complex[%d] %f %f\n", i, out_Complex[i].real(), out_Complex[i].imag());
	}
	for(i = 0 ; i < N ; i++){
		printf("out_cuda[%d] %f %f\n", i, out_cuda[i].real(), out_cuda[i].imag());
	}
	for(i = 0 ; i < N ; i++){
		printf("in_Complex[%d] %f %f\n", i, in_Complex[i].real(), in_Complex[i].imag());
	}
	for(i = 0 ; i < N ; i++){
		printf("out_cuda_inv[%d] %f %f\n", i, out_cuda_inv[i].real(), out_cuda_inv[i].imag());
	}

	fftwf_free(in); fftwf_free(out);

	getchar();
}
